Phase diagram and structural properties of a simple model for one-patch particles 
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We study the thermodynamic and structural properties of a simple, one-patch fluid model using 
the reference hypernetted-chain (RHNC) integral equation and specialized Monte Carlo simulations. 
Q . In this model, the interacting particles are hard spheres, each of which carries a single identical, 

CZ2 ■ arbitrarily-oriented, attractive circular patch on its surface; two spheres attract via a simple square- 

well potential only if the two patches on the spheres face each other within a specific angular range 
dictated by the size of the patch. For a ratio of attractive to repulsive surface of 0.8, we construct 
the RHNC fluid-fluid separation curve and compare with that obtained by Gibbs ensemble and 
grand canonical Monte Carlo simulations. We find that RHNC provides a quick and highly reliable 
estimate for the position of the fluid-fluid critical line. In addition, it gives a detailed (though 
g approximate) description of all structural properties and their dependence on patch size. 

o 

I. INTRODUCTION 

The role that anisotropic interactions play in the structure of a fluid system and on its equilibrium phase diagram 
qq | is the subject of intense current investigation; a renewed interest is reinvigorating previous studies that attempt to 
\Q i understand and model the structure and dynamics of associating liquids This renewed interest is stimulated 

t-H ' by two complementary recent developments. First, profound advances in chemical particle synthesis are significantly 
, augmenting the naturally available choices, offering the possibility of controlling surface patterns and functionalities in 
— i ' colloidal particles^ Command of the assembly properties of these new patchy colloidal particles will give fine control 
0^ , over the three-dimensional organization of materials as well as the combination of different materials over several 
■ length scales, promising to yield a spectrum of crystal polymorphs and assembled structures that is unprecedented 
| in colloid science. And second, growing evidence that to understand protein interactions in aqueous solution and to 
• ^ describe quantitatively their phase diagrams could require accounting for the presence of anisotropic protein-protein 
' interactions. Aeolotopioi^ interactions [from the Greek aiolos (variable) and topos (place)] arc now beginning to be 
^ . considered fundamental^^ for grasping the nature and mechanism of self-aggregation processes in protein systems, 
- - ' particularly in fibril and crystal formation. 

In the case of spherically symmetric potentials, a rather clear understanding has been reached of the topology of 
the phase diagram, of the structural properties of the liquid phase, and of the possible crystal structures. Much less 
understood are highly anisotropic models, where the local structure of such fluids is still rather poorly predicted. This 
is particularly true in those cases in which the anisotropy of the interaction or its patchiness produces first-nearest 
neighbor shells with low coordination number. Most of the evidence for the structure of these systems arises from 
numerical simulations , 11 i 12 ' 13 i 14 ' 15 i 16 which is only rarely accompanied by quantitative theoretical predictions. In this 
respect, primitive models of patchy particles^ — by condensing the interparticle repulsion into a hard-core potential 
while attractive interactions are modeled as simple square- well potentials — offer the possibility of establishing general 
trends and acting as test cases for checking novel theoretical approaches or extensions to the anisotropic case of familiar 
approximations whose validity has already been assessed for spherical potentials. 

One interesting primitive model is the Kern-Frenkel potential^ based upon an idea introduced by Chapman et 
al*£ in the framework of molecular fluids. In the simplest realization of this model, which we discuss in this article, 
the surface of the particle is partitioned into two parts, one with square-well and one with hard-sphere character. The 
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relative ratio between attractive and repulsive surfaces is measured as coverage x- Two particles attract if they are 
within a specific distance of each other and if their attractive surfaces are properly aligned with each other. 

Here we present an integral equation study of the phase diagram and of the structural properties of the one-patch 
Kern-Frenkel potential^ 7 - for a specific value of the coverage (x = 0.8) based on the reference hypernetted-chain 
(RHNC) approximation introduced in Ref. Ql^ for spherical potentials and later extended to molecular fluidsi 20 i 21 
This approach provides, besides thermodynamic information, an extremely detailed description of the structure of 
the liquid via a pair distribution function which now includes both radial and angular correlations, i.e., with a level 
of detail which can hardly be reached in experiments or in simulations. We compare the theoretical estimates of the 
gas-liquid coexistence and of appropriately-averaged radial distribution functions with "exact" numerical simulations 
of the model to assess the ability of the RHNC approach to describe such highly anisotropic interactions. 

As the size of the patch decreases, lower temperatures are required to follow the fluid-fluid coexistence line and 
integral equation convergence becomes progressively harder to achieve. Still, in the fluid region, it is possible to provide 
a detailed description of patch size dependence for the pair distribution function. We find a rich physical behavior 
which cannot be captured by rotationally averaged quantities, such as those computed by numerical simulations, nor 
by any effective temperature rescaling. 

Attractive interactions need not be restricted to short-range square- well potentials. In a very recent paper, Gogelein 
et al^ used a hard-sphere Yukawa potential, within a similar patchy model, to incorporate possible longer-range 
interactions due to the addition of salt, but they used second-order thermodynamic perturbation theory to derive the 
phase diagram, in contrast with the integral equation approach followed here. 

The paper is organized as follows. In Section [Til we introduce the one-patch Kern-Frenkel model used in this work, 
while in Section 1X111 we describe the main steps necessary for the solution of the RHNC integral equation of the model, 
stressing analogies and differences with similar approaches already exploited in the framework of molecular fluids; 
a few technical details are collected in the Appendices. Section llVI deals with the computation of thermodynamics 
and coexistence lines. Finally, Section [V] contains the RHNC results along with their comparison with Monte Carlo 
simulations and Section [VI] points out some conclusions and draws some perspectives. 



II. THE MODEL 



We study the Kern-Frenkel patchy hard sphere modelji 7 - patterned after an idea introduced by Chapman et al^ 
in the framework of molecular fluids: Two spherical particles attract via a short-range square-well potential only if 
two patches on the respective surfaces are properly aligned (see Fig. [JJ. 

In the case of a single patch per particle, the pair potential reads^ 



where 



and 



$(12) = <£(r 12 )*(ni,n2,?i 2 ), (1) 
oo, < r < a 

4>{ r ) = { <J <r < \<r (2) 
0, Act < r 

,j. » ~ \ f 1, if ni • fi2 > cos O and -n 2 • r 12 > cos O , f , 

*(m,n 2> n a ) = | 0) otherwise, (3) 

0o being the angular semi- amplitude of the patch. Here ni(wi) and n 2 (cj 2 ) are unit vectors giving the directions of 
the center of the patch in spheres 1 and 2, respectively, with lo\ = (9±,ipi) and w 2 = (0 2 ,ty5 2 ) their corresponding 
spherical angles in an arbitrarily oriented coordinate frame. Similarly, f i 2 (O) is the unit vector of the separation ri 2 
between the centers of mass of the two spheres and is defined by the spherical angle 0. As usual, we have denoted 
with fj and Aer the hard core diameter and the width of the well. 

Patchy models are the minimal version of models used to describe systems with highly directional attractions and 
have recently been applied to a variety of anisotropic cases, ranging from network and associated fluids j 4 ^ 23 ' 24 ^ 5 
colloids ) 14 ' 16 proteins ^ 17 ' 22 and gels. 15 ' 26 



III. SOLUTION OF THE ORNSTEIN-ZERNIKE EQUATION 



Our integral equation approach is based on the reference hypernetted-chain (RHNC) approximation introduced in 
Ref. [lj^ for spherical potentials and later extended to molecular fiuids i 20 i 21 In Appendix [XJ we briefly summarize 
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the necessary key expressions. Within this scheme, the closure equation takes on the assumed-known bridge function 
B Q (12) of a particular reference system to replace the actual unknown bridge function B(12) appearing in the exact 
closure (see Eq. (jA2|) ). Parameters of Bq{12) are treated as adjustable, a viewpoint initiated by Roscnfeld and 
Ashcroft(2i who postulate that the form of the bridge function B(r) for simple fluids constitutes the same univeral 
family of curves, irrespective of the pair potential, and that the known hard sphere version — with appropriate choice 
of the hard sphere diameter cto — can serve as a reference model for all of them. They demonstrate this hypothesis by 
explicit calculation^ selecting the optimum oo in each case so as to yield close agreement with published simulation 
results^ It is in fact possible to optimize the reference system parameters autonomously via a variational free energy 
principle 2 ^ that enhances internal thermodynamic consistency, in that Hiroikc's condition, 30 



dp\N J 8(3 \ p 

is recovered. This follows because, as in the original HNC equation, the internal energy U and pressure P are again 
obtained as thermodynamic derivatives of the same free energy functional. 

Anisotropic integral equations are far more demanding than their spherical counterparts from a computational 
point of view, but the basic procedure, involving expansions of the angular dependence in spherical harmonics, has 
been tested for numerous physical system s 31 ' 32 ' 33 with potentials having a continuous dependence on the angular 
orientations. The unusual discontinuous nature of the angular dependence appearing in Eq. ([3]), however, would 
seem to present a real challenge to such a spherical harmonics expansion, given especially that theoretically infinite 
expansions become necessarily finite in numerical work. In fact, there is no need to expand this potential at all; it 
may be used directly as-is within a specified angular grid. 

We briefly sketch now the main steps of the integral equation procedure. The iterative solution of the Ornstein- 
Zcrnike (OZ) equation plus an approximate closure equation typically requires a scries of back and forth Fourier 
transformations between functions in direct and momentum space. Besides the bridge function _B(12), this involves 
the direct correlation function c(12), the pair distribution and correlation functions <?(12) and h(l2) = g(12) — 1, 
respectively, and the auxiliary function 7(12) = h(12) — c(12) (see Appendix |A")) . 

Angular dependence introduces additional back and forth Clebsch-Gordan transformations between the coefficients 
of angular expansions in "axial" frames, with z = r in direct space or z = k in momentum space, and more general 
"space" frames with arbitrarily-oriented axes. (Note that here and in much of the following we have put for brevity 
f = fi2-) This is now a well-established procedure whose general scheme is summarized in Table I. We note that, 
because of the cylindrical symmetry of the angular dependence in Eq. ([3]) , the simpler version of the procedure for 
linear molecules^ - rather than the full anisotropic version 3 -^ can be exploited here. 

We start an iteration with the current values of the axial coefficients 7z 1 z 27Tt (r) in r space (see Appendix [A} and 
use the expansion in spherical harmonics, Eq. (|A3[) with X = 7, to construct 7(7", wi,^). The bridge function 
B(r,wi,u)2) is similarly constructed in RHNC approximation, while <5(r, u>\,u>2) is given. Then using the closure 
equation (|A2|) to obtain c(r, Wi,^) and the expansion inverse Eq. (|A4[) (with X = c) to compute its axial coefficients 
C; 1 ; 2m (r), the algorithm proceeds by a Clebsch-Gordan transformation to obtain the corresponding space coefficients 
c(r; hfol) from Eq. (|A6[) . This is followed by a Hankel transform, Eq. (|A8|) . to momentum space and a backward 
Clebsch-Gordan transformation, Eq. (|A7[) . to return to an axial frame in k space, yielding ci 1 i 27n (k). The OZ equation 
(1A5[) in momentum space can now be invoked to get ^yi^^m 

(k). This is then followed by a forward Clebsch-Gordan 
transformation, Eq. (|A6[) . and an inverse Hankel transform, Eq. (|A9[) , to provide j(r;l±l2l). A final backward 
Clebsch-Gordan transformation, Eq. (|A7[) . brings a new estimate of the original coefficients 7z 1 z 2m (r). 
This cycle is repeated until self-consistency between input and output 7i 1 z 2 m(»") is achieved. 



IV. THERMODYNAMICS 



A. Internal energy and virial pressure 

The internal energy U of the system can be computed starting from the standard expressior 

N N 
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u = ( E E $ (y) ) = \ pN I dri2 {9 (12) $ (12)) — ■ (5) 



i i— 1 j>i 



Here p is the number density and (. . .) u = J du . . . indicates an average over the spherical angle u> = (9, tp). 

This expression can be written in a simpler form by exploiting the factorization in Eq. ([1]): 

^ = lp J dr (3<t> (r) G (r) = -2npf3e J dr r 2 G (r) , (6) 
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where (3 = 1/ksT is the inverse Kelvin temperature (with ks Boltzmann's constant) and 

1 



G(r) = (g(r, ui, w 2 ) * {uji, w 2 )) Wl w 2 = 



(4^ 



du>iduj2 g (r, Ui, u>2) * 0J2) 



(7) 



The angular integrals in Eq. |(7J) are evaluated using Gaussian quadratures^ while the trapezoidal rule is used for 
Eq. ©. (The average number of neighbors in the well is found as z = —2U/Ne.) 
The virial pressure P can also be obtained from a standard expression^ 



^ I N N q \ 1/" / \ 

P = ^ BT -3y\EE^^T $ W)y =pfcsT--p 2 y dr 12 ^(12)r 12 — $(12)^. 



(8) 



Using familiar manipulations, the definition of the continuous cavity function y(l2) = g(12)e^^ 12 \ and the result 
d 



3 -/30(ri 2 )*(wi,w 2 ) 



dri2 

Eq. ([5]) can be cast in the form 



3 /3e*(wi,w 2 ) 



1 



S (ri2 - Aa) 



7" = i + ^ 3 < 



- 1 



(9) 



(10) 



which again can be reduced to Gaussian quadratures. 

Finally, the isothermal compressibility kt = p~ 1 (dp/dP)x can be obtained from^ 3 . 



pk B Tti T = 1 + p J dr [5000 (r) - 1] 



(11) 



B. RHNC free energy 

It has been shownSi that in RHNC approximation, just as with its parent HNC equation^ the excess free energy 
per molecule can be calculated in a self-contained fashion as 

(3F CX [3F 1 , 0F 2 , 0F 3 



N 



N 



N N 



where 



N 



~p [ dr 12 /hi 2 (12) + h(12)-g(12)ln [g (12) e ^ 12 > 



= -tI 7^3E{ lnDet f I +(- i r^ TO (fc)i -(-ir P Tr\h m (k)}}, 

N 2pJ (2wY _ L L J 



(2^) J 



(12) 

(13) 
(14) 

(15) 



In Eq. (|14p . h m (fc) is a Hermitian matrix with elements /i; 1 i 2m (fc), /i,Z 2 > m, and I is the unit matrix. The last 
equation, for F3, directly expresses the RHNC approximation. Here F3 is the reference system contribution, computed 
from the known free energy F e ° x of the reference system as F® = F e ° x — Fj 3 — F^, with F® and F® calculated as above 
but with reference system quantities. 

The bridge function i?o(12) appearing in (|15|) is taken from a reference system and is the key approximation in 
the RHNC scheme, since it replaces the unknown bridge function -6(12) in the general closure equation (|A2[) . At 
present only the hard sphere (HS) model is sufficiently well known to play the role needed here, and so we put 
Bo (12) = -E>hs(?"i 2 ; Co). The optimum hard sphere diameter 00 is selected according to a variational free energy 
minimization that yields the condition 2 ^ 



p J dr [5000 (r) - gas (r; cr ] 



(9-Bhs ( r ;°o) 
dcr 



= 0. 



(16) 



Note that cr = corresponds to an HNC closure while cr = 1 corresponds to the true hard-sphere diameter of the 
potential. In actual numerical computations, the value of a$ clearly will depend upon the considered state points and 
coverage. This dependence of ctq in our results is further discussed in Section [VB1 
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C. Chemical potential 



The chemical potential p is conveniently obtained within the RHNC scheme from the relation 

^= W + T , (17) 

which does not introduce any additional approximation into the calculation. The ideal quantities are PF-^/N = 
ln( j oA 3 ) — 1, pPia/p = 1, and /3/i;d = ln^pA 3 ), where A is the de Broglie wavelength. 

We stress the fact that, within the RHNC scheme, the chemical potential directly follows from the free energy 
and pressure, themselves computed directly from the integral equation solution for the pair functions, unlike other 
closures where additional approximations are required. This is a particularly important feature when numerical results 
stemming from different closures are compared. 



V. NUMERICAL RESULTS 



A. Monte Carlo simulations 



We have studied the equilibrium properties of the one-patch model with several standard Monte Carlo methods. We 
perform Gibbs ensemble simulations 3 ^ to evaluate coexistence in the region where the gas-liquid free-energy barrier 
is sufficiently high to avoid crossing between the two phases. Here we use a system of 1200 particles which partition 
themselves into two boxes whose total volume is 4300er 3 , corresponding to an average density of pa 3 = 0.27. At the 
lowest T this corresponds to roughly 1050 particles in the liquid box and 150 particles in the gas box (of side ~ 13er). 
On average, the code attempts one volume change every five particle-swap moves and 500 displacement moves. Each 
displacement move is composed of a simultaneous random translation of the particle center (uniformly distributed 
between ±0.05ct ) and a rotation (with an angle uniformly distributed between ±0.1 radians) around a random axis. To 
calculate the location of the gas-liquid critical point, we perform grand canonical Monte Carlo (GCMC) simulations^ 
complemented with histogram reweighting techniques to match the distribution of the order parameter p — se with the 
known functional dependence expected at the Ising universality class critical points Here e is the potential energy 
density, p the number density, and s the mixing field parameter. We have studied systems of size 7a up to 10cr to 
estimate the size dependence of the critical point. We did not perform a finite-size study since the changes in critical 
parameters with the size are not significant. In GCMC we attempt, on average, one insertion/deletion step every 500 
displacement steps. We have also performed a set of GCMC simulations for different choices of T and p to evaluate 
the corresponding density. Finally, we performed standard MC simulations of a system of 1000 particles in the NVT 
ensemble to estimate structural properties for comparison with integral equations results. 



B. Integral equations 



We solve the RHNC equations numerically on r and k grids of N r = 2048 points, with intervals Ar = 0.01a and 
Afc = Tr/(N r Ar), using a standard Picard iteration method^ The square-well width is set at A = 1.5 because of 
results available for the isotropic square well.— For the hard-sphere bridge function, we use the Verlet-Weis-Hcnderson- 
Grundke paramctrization . 39 ' 40 Reduced temperature T* = ksT/e and density p* = pa 3 arc used throughout. 

The angles 6j, j — 1, 2, . . . , n, that appear in the Gauss-Legendre quadratures over polar angle 9 are determined by 
the n zeroes of the Legendre polynomial P^cos^) of order n; it is not possible in general to ensure that one of these 
will coincide with the patch parameter 8q. To obtain a "best fit," we have chosen n in the following ad hoc fashion. 
The total coverage x (fraction of surface covered by the attractive patch) can be obtained from Eq. as (cf. Ref. 

m 

X = (^(ftx.n^ria^^aiii 2 ^. (18) 

By also numerically evaluating the angular average in this expression using Gauss-Legendre quadrature of order n 
and varying n (typically from 30 to 40), we find that number of Gaussian points n that minimizes the error in the 
numerical computation of x, the exact value being known. Then for a fixed value of x, all Gaussian quadratures in 
a calculation are evaluated with the same number n of points. For the case x = 0-8 studied in most detail in this 
article, we found n = 31 to be the optimal value, yielding a 0.25% error in Eq. (|18p . For lower x values, the value 
n = 30 was used instead. 
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Expansions in spherical harmonics are terminated after l\ = I2 = m = 4, resulting in 35 distinct coefficients in a 
pair function expansion; we have explicitly checked that this suffices for all cases considered here. Hankel transforms 
have been obtained using a step- up/step-down procedure as detailed in Ref. HO. In all calculations, the covergence 
criterion was a root-mean-square difference between successive iterates of less than 10~ 5 . At 80% coverage, typical 
values for the variational parameter oq range from approximatively 0.77 for reduced densities of order I0 -2 to 1.02 
in the liquid branch (reduced density 0.75). These latter values decrease for lower coverages to reach the HS limit f .0 
when x — > 0, as expected. 

C. Phase coexistence 

For the computation of the phase diagram we have here considered a representative value of coverage \ = 0-8 
Figure [2] depicts a comparison of the results for the reduced density p* as a function of the activity e 13 ^. Here p is the 
total chemical potential, so that p* = e 13 ^ is the ideal gas limit, where we have replaced the density- independent de 
Broglie wavelength A with the hard sphere a for convenience. A comparison with GEMC results shows that RHNC 
reproduces the correct behavior both qualitatively and quantitatively for both the low (gas) and high (liquid) densities 
branches. This is a rather stringent test on the performance of RHNC since the chemical potential is usually a very 
sensitive quantity to the approximation used. For the liquid branch, the numerical values stemming from the two 
procedures are also reported in Table HT1 

Armed with this result we turn next to the computation of the phase diagram. There exist very few studies of 
the fluid-fluid coexistence line using integral equations even in the framework of spherical, fully isotropic, potentials. 
One possible scheme to identify the densities of the two coexisting phases was outlined by Schlijpcr et a/4i for a 
Lennard- Jones fluid. A similar method was used in Ref. [38| for a square-well fluid. A more elaborate procedure has 
been devised in Refs^ 2 - and^. Our analysis closely follows the spherical counterpart given in Ref. [H in order to make 
a close connection with the fully attractive case. 

For a fixed temperature below the critical one, the gas and liquid branches are scanned along the density axis until 
the numerical program is unable to produce a converged solution. This might be due either to the disappearance 
of a solution of the integral equation itself, which is a well-known shortcoming^ of the hypernetted-chain closure 
and its descendants, or to the approach of spinodal instability. Coexistence density values p g and pi for the chosen 
temperature are then obtained by the numerical solution of the system 

P{pt) = P(P S ), (19) 

fi(jn) = n(p g ). (20) 

It might (and it docs) happen that approaching the spinodal line, it is not possible to obtain a solution of the 
integral equation before Eq. (jf 9[) can be satisfied. In such cases, a suitable small extrapolation compensates for this 
limitation of the approximate closure; details of this extrapolation process can be found in Ref. \3§L Figure [3] depicts 
the phase diagram so obtained compared with GEMC results, where extrapolated points have been highlighted. Wc 
note that RHNC results follow remarkably well GEMC data along both branches. For RHNC, the numerical values of 
the (reduced) densities and temperature along with the common values of chemical potentials and pressure, as given 
by Eq. (fT9]) . are listed in Tabic HTT1 The small systematic deviation from the GEMC results in the high density branch 
is roughly of the same order of magnitude as that found in the SW case^ thus strongly suggesting that angular 
dependence has not worsened the results. This discrepancy would obviously be reduced by a better description of the 
anisotropic bridge function, about which essentially nothing is known. Minor deviations in the low density branch, on 
the other hand, arc mostly due to numerical convergence problems, especially for temperatures close to the critical 
region, as indicated by the extrapolation marks. 

As the coverage \ decreases, we find that the critical temperatures and densities shift towards progressively smaller 
values. This happens also for higher number of patches and in general it cannot be accounted for by a simple rescaling 
of the temperature^ We remark, however, that the one-patch case is very peculiar in this respect for two reasons. 
First, because it is the only case (along with its two-patch counterpart) clearly interpolating between the full isotropic 
HS (x = 0) and SW (x = 1) limits. For 3 or more patches on each sphere, a well-defined coverage limit (smaller than 
100%) is found. Second, because of new peculiar phenomena arising for low coverages, notably for \ < 0.5. This and 
other points will the the subject of a future publication. 

D. Pair correlation functions and structure factor for \ = 0.8 coverage 

Given the results obtained for the phase diagram, it is interesting to consider structural information from an 
analysis of correlation functions and structure factors. Unlike with their isotropic counterparts, information about 



7 



pair distribution functions g(12) with angular dependence is obtained differently in integral equation calculations 
and Monte Carlo simulations. Integral equations provide a complete, although approximate, description of the pair 
distribution function g(12), which in the "canonical" axial frame is a function of four variables: ri2, 0i, 02, and 
012 = $2 — 4>i- This level of detail is a prohibitive computational task for simulations, which instead produce a more 
accessible average of some sort. In the present simulations, we calculate an average over all possible orientations of 
the vector joining the centers of mass of the two spheres, fi2(0), while holding fixed the orientations of the individual 
vectors ni(wi) and n2(<x>2), as well as the magnitude ri2- This produces a function of just two variables: 7*12 and 

012 =#2—01- 

Detailed results for specific orientations from RHNC, without averaging, can be seen in Fig. [U where we plot 
g(12) = g(r, u>i,u>2) as a function of r/a for two state points lying in the liquid-like region (p* = 0.68) with two 
different temperatures: T* = 1.0 above the critical point and T* = 0.65 in the subcritical region. The three different 
curves in each plot refer to three different values of the product 111-112; these are ni • ng = 1 for a head-to-tail (HT) 
configuration, for a crossed (X) configuration, and -1 for a head-to-head (HH) configuration. 

Coverage has again been set at 80% and the width of the well is as always A = 1.5. The jump occurring at r/a = 1.5 
in the HH and X configurations closely resembles the jump occurring in the isotropic counterpart^ indicating a clear 
bonding due to the attractive part of the potential. This does not occur for the HT configuration, where the patches 
are on opposite sides (see Figj2]). This feature is common to all coverages smaller than 100% (i.e., a square- well 
potential) and it is related to the particular angular dependence of the Kern-Frenkel model, as given in Eq.Q, 
coupled with the fact that g(12) is here computed in the axial frame where fii is along the direction f 12 . Then 
Ai • ii2 = 1 and — A2 • V12 = —1 in the HT configuration thus yielding no coupling unless 9q = 7r corresponding to 
the isotropic square- well potential. Results corresponding to the isotropic square- well potential (100% coverage) and 
hard-sphere potential (0% coverage) are also depicted for contrast. 

A few features of the above behavior are worth noting. First the hard core contact values (r = a + ) in the HH and 
HT configurations are virtually identical at both temperatures T* = 1.0 and T* = 0.65, whereas the X configuration 
is significantly different, this difference being greater for the lower T* = 0.65 temperature. 

As temperature decreases from T* = 1.0 to T* = 0.65, the SW contact value decreases, the X contact value 
also decreases while the corresponding HH and HT ones are essentially unchanged. This directionality ordering is 
reversed at the other contact point inside the well (r = Act~) where the SW value is larger than both the HH and 
X corresponding ones, and this difference increases with increasing temperature. The HT configuration has a much 
lower value as there is no discontinuity in this case, as remarked. Inspection of the HT contact values is sufficient 
to convince that the head-to-tail configurations display a significant temperature dependence in spite of their hard- 
spherc-likc apparent behaviors. 

The above results can be understood by contrasting these cases with attractive interactions with the corresponding 
HS ones (independent of temperature) . Here the r = a + contact value is larger than all other cases (including the SW 
one) where the other r = Act - extreme has a much lower value compared to SW, HH and HT configurations. This is 
because attractive interactions decrease the probability of sphere-touching occurrence while it increases the probability 
of an indirect bonding within the well region a < r < Act mediated by the presence of intermediate particles. For 
lower coverage, the overall attractive force is reduced and this effect is less marked and, more importantly, it depends 
upon the directionality of the bond. 

A corresponding low-density (p* = 0.1, T* = 1.0) state is also reported in Fig. |4]in order to account for a gas-like 
behavior. Not surprisingly, in this case both HH and X configurations behave similarly, with HH configuration always 
lying above the corresponding X one, since for such a low density the low probability of finding nearest-neighbor 
spheres is mostly independent of the orientation (the fully isotropic results is not reported in Fig[4]as this state point 
lies inside the critical line; see Ref. [HI). 

Snapshots of configurations for the above state points are displayed in Fig. [5J The liquid-like structure is clearly 
visible in the first two configurations with p* = 0.68, without any significant morphological difference between the 
supercritical T* = 1.0 and the subcritical T* — 0.65 thermodynamic states. The RHNC provides also a prediction 
for the average number of particles in the well z (that is, with center-to-center distance within the range [ct, Act]) by 
integrating the center-to-center radial distribution function (g(12)) Wl:U2i n within the well. In the above cases we find 
approximatively 7.5 and 8.3 for T* = 1.0 and T* = 0.65, respectively. 

Conversely, in the vapor- like state p* =0.1, T* = 1.0, a much lower average number of particles, 1.6, is found in 
the well, as expected. This figures have been obtained upon integration of the g(12) over the radial range a < r < Xr 
averaged over all possible relative orientations. 

Since the integral equation results provide a complete, albeit approximate, picture it is possible to calculate from 
them the simulation average of <?(12) over all orientations of f 12(^1), while holding ni, n.2, and r\2 fixed. This is 
detailed in Appendix [Bl The result is a quantity we denote g(r, cos0i2), which can be directly compared with the 
Monte Carlo results; this is done in Fig. [51 An overall agreement between MC and RHNC results can be seen, although 
the RHNC performance clearly worsens as temperature decreases. It is worth emphasizing that, upon averaging, even 
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an HT configuration shows a jump at r/er = 1.5, an indication of bonding in some orientation of fi2(0). 

Additional information about structural properties can be obtained from structure factors. In the presence of angular 
dependence, the quantity playing the role of the usual center-to-center structure factor is the Sooo(k) coefficient (see 
Appendix [A| of S(12) = 1 + ph(12), which is plotted in Fig. [7] as a function of ka for the same state points as 
above. This spherical component of the structure is rather well predicted by the theory. Viewing the temperature 
dependence of Sooo(k), one notices that on cooling the second peak splits into two distinct peaks. Such a splitting, less 
evident but present also in the fourth peak, can be traced back to interference between the asymptotic structure factor 
oscillations due to the two discontinuities of the pair correlation functions at a and Act. By decreasing temperature, 
the contribution of the pair correlation function close to Act becomes comparable with that from the region close to 
the hard sphere contact. 

E. Pair correlation functions and structure factor at different coverages 

As coverage decreases, it becomes harder and harder for integral equations to follow the critical line. This is because 
a lower coverage corresponds to a lower overall attractive interaction. As such, the critical temperatures and densities 
progressively decrease and integral equations approximations and algorithms become less reliable. 

Yet we can obtain interesting guidelines of the general trend by considering structural properties at different 
coverages. We then select a specific state point p* = 0.68 and T* = 1.0, which is above the critical temperatures 
for all coverages but still displays significant structural behavior. Table IIVI reports the RHNC results for the main 
reduced thermodynamic quantities previously defined: the internal energy PU/N per molecule, the excess free energy 
per molecule (3F eK /N, the chemical potential [3p and the compressibility factor (3Pj p. In addition we have also 
calculated the quantity /3(dP/dp), whose approach to zero signals the vicinity of a spinodal line. 

Figure [S] again depicts g(12) in the HH (top panel), X (central panel), and HT (bottom panel) configurations as a 
function of the coverage starting from a fully square- well potential (x = 1-0) to a hard-sphere potential (x = 0.0). The 
non-monotonic dependence on coverage at the SW discontinuities r = a + and r = Act - is particularly noteworthy and 
has been singled out in Table fVl For brevity, we have put there g a {r) = g{r,u>x,u)2) for r = ct + and r = Act - , where 
a stands for HH, X, or HT. Contact values within the well follow a rather erratic pattern as detailed in Table Ivl and 
Fig. [SI The contact value <7hh(c + ) first increases as \ decreases from 1 to 0.5 and then oscillates until a final sharp 
decrease to the hard-sphere limit at \ = 0-0. The particular high value of \ — 0-1 ( sce Table |V|) is a clear indication 
of a chaining and micelles formation, as one expects for such a low coverage, but similar - albeit lower — values are 
found for \ = 0.2 and x — 0-4, always below the 50% limit. An even more irregular behavior can be observed at the 
other extreme of the well, 5hh(Act - ). 

This rich complexity cannot be explained by simple geometrical arguments hinging upon rotationally-averaged 
quantities. This can be quickly inferred by inspection of Tabic fVTI where we report the average correlation function 
g{r, cosO) defined as in the previous subsection, at the same interior well extremes ct + and Act - . Figure [5] also depicts 
the coverage dependence of g(r,—l), g(r,0), and g(r,+l), respectively HH, X, and HT. Clearly, previous irregular 
behavior in the g a (r), a = {HH,X,HT}, have been smoothed out and a much more regular trend is observed. For 
instance (see second column of Table IVT]) . <?(ct + ,— 1) gradually increases as x decreases from 1.0, peaks at 0.3, and 
drops to the hard-sphere limit. This non-monotonic behavior cannot be explained by simple heuristic arguments 
which would suggest a monotonic trend. As x decreases, the whole coexistence curve is lowered on the temperature 
axis, effectively corresponding to an overall decrease of the attractive forces or, equivalently, to an increase of the 
temperature scale. This is consistent with a decrease of the contact value of the full coverage (square-well) contact 
value (sec Figf?]top and central panel). Similar reasoning explains also the opposite trend followed by <?(Act - , —1) 
(third column in Table PVTl again contrasted with the square- well Act - behavior in Fig. 2] top and central panel, which 
increases as temperature decreases). 

The perpendicular X configuration follows a similar pattern (see Fig. [5] and fourth and fifth columns in Tables IVl 
and IVI[) but has a sharp change in the g(l2) behavior upon passing from x > 0.5 to x < 0.5 (see Table IVf indicating 
that 50% coverage plays a special role, as one could expect from the the outset. Again these nuances are lost upon 
rotational averaging (see Table [*VT|) . 

Next, we note that HT <?(12) configuration does not display any discontinuity at the well edge Act - (sec bottom 
panel of Fig. |5J) as already remarked. Nevertheless, there is a significant coverage dependence interpolating between 
square-well and hard-sphere behavior again in a non-monotonic fashion which is smoothed out in corresponding 
rotational averaged quantity (see Fig[5]and last two columns in Tables IVl and IVTf . 

As a final point, it is instructive to consider the coverage dependence of the radial structure factor Sboo(fe)- This 
is reported in Fig. [TOl As coverage decreases 5ooo(0) shows a monotonic decrease, in agreement with the fact 
that the critical region is shifting to lower temperatures. The first peak is slightly shifted to higher ka values and 
decreases as coverage decreases. The second peak displays an opposite variation with the coverage, accompanied by 
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an increasing width. Such behavior can be understood as due to the increasing contribution of the discontinuity at 
Act with increasing x (see discussion of the interference of the asymptotic oscillations of the structure factors at the 
end of Section V D). 

VI. CONCLUSIONS AND OUTLOOK 

We have studied the Kern-Frenkel model for hard spheres with a single attractive patch using integral equation 
techniques within the RHNC approximation and have compared results with various specialized Monte Carlo simula- 
tions. This potential constitutes a serious challenge to integral equation approaches for several reasons. First, because 
of the complexity in the underlying theory, which far exceeds that of its symmetrical counterpart. This is mirrored 
in the increased complexity of the corresponding numerical algorithms (see Appendix [A}. Second, the discontinuous 
nature of the patchy anisotropy poses in general a serious problem for theories based on expansions in spherical 
harmonics, as an infinite number of terms is in principle required to describe the discontinuity. While the first point 
is common to all molecular fluids, this second point is specific to patchy potentials. Remarkably, in the Kern-Frenkel 
potential case, it is possible to cope with this in a controlled manner by simply adapting the angular grid as detailed 
in the text. Finally, care is required when comparing pair correlation functions with Monte Carlo simulations, as 
different types of averages may be employed in simulations. Yet, since integral equations provide essentially complete 
descriptions within a given approximation defined by the closure, the pair correlation averages typically computed in 
MC simulations can be conveniently extracted from their results (see Appendix [B| . 

The RHNC integral equation has already been used in the past in the framework of associating molecular fluidsi 31 i 32 
It provides a robust and reliable description for a number of anisotropic potentials such as multipolar dependence. A 
major advantage of this particular closure, rendering RHNC particularly suited for computation of phase diagrams, is 
that the chemical potential calculation requires no additional approximation, thus enhancing internal thermodynamic 
consistency. This contrasts strongly with all other standard integral equation approaches (aside from the parent 
HNC), where pressure and chemical potential involve different approximations. Given a plausible representation of 
the bridge function, everything else follows, thus providing much better control over the procedures used. 

We have given a full description of the phase diagram and structural properties of the one-patch model with 
a specific coverage (80%). Notwithstanding a small systematic deviation in the detailed values, molecular RHNC 
correctly reproduces the general quantitative features of the phase diagram with an accuracy comparable to the 
analogous calculation of its spherical counterpart - a simple square-well potential. 

The pair correlation functions have been calculated for representative state points in the phase diagrams for three 
specific relative orientations of the patches on the pair spheres: when the two patches are aligned along the same 
direction (head-to-tail configurations), when they form a right angle (crossed configuration), and when they face each 
other (head-to-head configuration). The jumps characteristic of the square-well potential, the contact values, and the 
depth and position of the oscillations, all can be put in direct and natural correlation with the peculiarities of the 
patch dependence. Likewise, the k = value of the structure factor 5*000 (k) an d the position and morphology of the 
principal peaks can be directly related to the aggregation properties of the fluid and an appropriate angular average 
of the pair correlation function is amenable to a direct comparison with Monte Carlo simulations. We have shown 
that, even in the fluid-like regime, RHNC provides a rather precise description. 

Finally, we have provided a detailed study of the pair correlation functions and structure factor for different values 
of the coverage (or, equivalently, for different size of the patch). We find a flourish of different behaviors, with a 
non-monotonic dependence on the coverage, which are inaccessible to numerical simulations based on rotationally 
averaged quantities and could help to interpret the low-coverage fate of the coexistence line. 

The present RHNC study has two main limitations. First we observe a significant, albeit small, quantitative 
difference with numerical simulations in the prediction of the critical line location. This is due on the one hand to 
intrinsic drawbacks of the RHNC approximations which are present also in the fully isotropic spherical case^& and, on 
the other hand, to the extremely anisotropic nature of the single-patch potential, which emphasizes the shortcomings 
of the isotropic approximation used for the bridge function in the present work. Secondly we have restricted to a single 
representative coverage (80%) for the coexistence line location. This is because of the numerical instabilities arising 
for lower coverages in view of the much lower temperatures required. Both those shortcomings will be addressed in a 
future publication. 

A number of possible perspectives can be envisaged from the results of the present paper. As the coverage decreases, 
the binodal curve shrinks and moves to the left, thus rendering integral equation algorithms less and less stable, as 
remarked. This is nevertheless an extremely interesting region to study with integral equations, due to the subtle 
interplay between condensation and polymerization which could give rise to a very rich phenomenology 45 ' 46 ! 47 In 
addition, a comparison with the case of two patches, set on opposite sides of a sphere, at the same global coverage 
could illuminate the critical-point dependence on both valence and size distribution, a question of experimental 
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relevance. Preliminary results for the two-patches case clearly show that indeed RHNC is able to closely follow 
numerical simulation predictions for the critical line at different coverages. A detailed analysis of this and other 
related points will be reported elsewhere. 
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APPENDIX A: KEY EXPRESSIONS FOR THE INTEGRAL EQUATION PROCEDURE 

For the sake of compactness and future reference, we provide below the main expressions necessary for the numerical 
solution of OZ-type integral equations with orientational dependence. The resulting iteration scheme is laid out in 
Table I. 

1. Ornstein-Zernike equation and closure 

The OZ equation in terms of 7(12) = h(12) - c(12) is 

7(12) = £ J dr 3 ^ 3 [7(13) + c(13)] C (32), (Al) 

while the general form of a closure is 

c(12) = exp[-/?$(12) +7(12) +£(12)]- 1 - 7 (12) . (A2) 

The bridge function B (12) is known only as an infinite power series in density whose coefficients cannot be readily 
calculated^ All practical closures approximate B (12) in some way. 

2. Expansions in spherical harmonics 

The expansion in spherical harmonics Yi m (u>) of any correlation function A(12) in its axial frame (z = f) is given 

by 

X(12) = X(r,« 1 ,w 2 )=47r ]T X hhm (r) Y hm Y hfh {u 2 ) , (A3) 



where m = —m. and its inverse reads 



1 



X hhm (r) = — I <k> x du 2 X (r, u x , w 2 ) Y* m (uj Y*^ (w 2 ) . (A4) 

The same expressions, with k replacing r, apply to the Fourier transform A(12) in momentum space, where the 
axial frame is now defined by z = k. The coefficients in the expansions satisfy the symmetries X^^m = X^i.^ and 
X hhm = (-l) h+h X hhm . 

The angular integrations in the expansion inverse (|A4[) are evaluated using Gaussian quadratures. Note that to 
compute Q 1 i 2TO (r) from the closure equation (|A2[) using (|A4[) . we must first construct 7(12) and (an approximate) 
B(12) from their axial expansion coefficients using (|A3|) : the potential $ (12) is already in angular representation. 

The OZ equation factors in momentum space; expressed in terms of the axial expansion coefficients of the trans- 
formed pair functions, we get the simple matrix form 

00 

lhl 2 m(k) = (-l) m [ji lhm {k) + d hhm {k)}di 3 i 2 , n {k) 7 (A5) 

which can be solved for the 7i 1 ; 2m {k) using standard matrix operations. (A carryover from the Hankel and Clebsch- 
Gordan transforms defined below is that transform coefficients with li + Z 2 even are real while for li + 1 2 odd they are 
imaginary. Matrices such as c m (k) with elements Q 1 i 2m (fc), li,l 2 > m, are Hermitian.) 
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3. Clebsch-Gordan transforms 



For the function X(12) in direct space, the Clebsch-Gordan transform from axial coefficients X; 1 i 2rra (r) when z = f 
to space coefficients X (r; l\l 2 l) is 

4tt ^ /2 - 



x ^ khk = lu + i) E c W;mmo)x W2TO (r), ( A6 ) 

where the C (hl 2 l; mfhO) are Clebsch-Gordan coefficients; the inverse transform is given by 

/2/ 4- 1 \ 

*W a m(r) = ^Cfl^mmO) ^ X(r;hl 2 l). (A7) 



The same Clebsch-Gordan transforms, with A: replacing r, apply to the Fourier transform X(12) in momentum space, 
where the axial frame is now defined by z = k. Symmetry requires that l\ + l 2 + 1 be an even integer. 



4. Hankel transforms 



The Fourier transform in an arbitrary space frame of a function X(12) into X(12) leads to the Hankel transforms 
of the coefficients, 

/>OC 

X(k;hhl) = 47ri' / dr r 2 X (r;hl 2 l) ji (kr) , (A8) 
Jo 

with the inverse transform reading 

X(r;hl 2 l) = -^r / dkk 2 X(k;hhl)jl{kr), (A9) 
^ 1 Jo 

where 7z(a;) is the spherical Bcssell function of order /. By using "raising" and "lowering" operations on the 
integrands^ these are finally evaluated with jo(kr) = smkr/kr kernels, for I even, and j~i{kr) = coskr/kr ker- 
nels, for I odd. Fast Fourier Transforms are programmed for both instances. (Note that Hankel transforms X (k; l\l 2 l) 
are imaginary for I odd while for I even they are real.) 



APPENDIX B: COMPARING INTEGRAL EQUATION PAIR CORRELATION FUNCTIONS WITH 

NUMERICAL SIMULATION RESULTS 

The function X(12) 1 expanded in its axial frame where z = f, is given in Eq. (|A3[) . If we rotate the frame of 
reference into an arbitrary space frame, its spherical harmonic expansion becomes^ 

X(U) = Air X ( r -Mhl) C (y 2 i; m 1 m 2 m 1 + m 2 ) 7 iirai }^ 2m2 ( W2 ) y ; ; mi+m2 (O) , (Bl) 

(1,(2,2 mi, ma 

where the X (r; l\l 2 l) are related to the original X; 1 / 2m (r) through the Clebsch-Gordan transform Eq. (|A6|) . In order 
to compare correlation functions such as g(12) with the results of our Monte Carlo simulations, we need to average 
g(12) written as in Eq. (|B1|) over all possible orientations of f (fl), keeping ni(wi), h 2 (uj 2 ), and r fixed. By symmetry, 
the result can depend only on r and cos #12 = fti • £12. We define thus 

g(r,cos6 12 ) = i- f dQ g{12). (B2) 

We find first that 

^ / <m Y£ mi+ m 2 (n) = j^yj- 2 J dn r„ (0) Y i:mi+m2 (0) = SmS £ffi\ (B3) 
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so that Eq. (|B2j) becomes, with (|BT 



Then using^ 



3 (r, cos 6*i2 2?n (^2) • (B4) 

li ,12 ,m 



C(hl 2 0;mmO) = l^ 1)1/2 ftifa (B5) 
33 



and the addition theorem for spherical harmonics,— one finds finally that 



°° /?Z + 1\ 1/2 

ff(r,cos0 12 ) = £(-1) I "tH 9 (r\U0) Pi (cos 9 12 ), (B6) 



where Pi(cos8) is the usual Legendre polynomial of order I. 
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, « Hankcl transform Eq. JA8t ~ / ? j\ Inverse CG transform Eq. jA7t „ , . 

c{r;hl2l) > c(fc;/i/ 2 () ► ci^i 2m (k) 

JcG transform Eq. fACl |oZ equation Eq. fAsl 

J Expansion inverse Eq. iA4\ j,*^^ trans ^ orm ■^ c ^■ 

c(r, wijWa) l{k\lihl) 

J Closure Eq. \A2\ ^Inverse Hankcl transform Eq. JA9t 

7(r,£Ji,a;2) 7(r; J1/2Z) 

^Expansion Eq. j A3| I Inverse CG transform Eq. ]A7| 

[ihhm (r)] old < Iterate < [ihh™ Ml™ 

TABLE I: Schematic flow-chart for the solution of the OZ equation for angle-dependent potentials. 
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exp(P(jt) (GEMC) 


p* (GEMC) 


exp(/3p) (RHNC) 


p* (RHNC) 


0.0105 


0.6044 


0.0119 


0.65 


0.0120 


0.6287 


0.0165 


0.68 


0.0150 


0.6555 


0.0219 


0.70 


0.0180 


0.6752 


0.0380 


0.73 


0.0210 


0.6845 


0.0596 


0.75 


0.0240 


0.6934 


0.0997 


0.77 






0.2408 


0.88 



TABLE II: Tabulated values for exp((3p) and p* corresponding to the liquid branch of Figure [5] for both GEMC and RHNC 
data. Note that these values are results of different calculations in the GEMC and RHNC cases and the two abscissas are not 
the same. Here and below expected errors are in the last digits. 
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T* 


P* 3 


Pi 


Pn 


PP 


0.50 


0.0016 


0.7465 


-6.4789 


0.0016 


0.55 


0.0040 


0.7173 


-5.6342 


0.0038 


0.60 


0.0086 


0.6870 


-4.9607 


0.0077 


0.65 


0.0171 


0.6518 


-4.4147 


0.0141 


0.70 


0.0306 


0.6098 


-3.9668 


0.0228 


0.75 


0.0414 


0.5436 


-3.5965 


0.0340 



TABLE III: Numerical values of the RHNC coexistence line of Figure [3] Here p* and p* axe the reduced densities of the gas 
and liquid branches. T* is the corresponding coexistence reduced temperature, and Pp and PP are the common values of the 
chemical potential and pressure. 
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Y 




u 

JVe 


/3-Fcx 
JV 




p 






£0 
a 


1. 





-5.32 


-2.55 


-3.58 


0.35 


8.62 


1 


.031 





9 


-4.52 


-1.70 


-1.82 


1.27 


9.87 


1 


.025 





.8 


-3.76 


-0.92 


-0.24 


2.07 


10.86 


1 


.018 





7 


-3.02 


-0.18 


1.24 


2.81 


11.98 


1 


.014 





6 


-2.37 


0.45 


2.49 


3.43 


12.88 


1 


.010 





5 


-1.77 


1.02 


3.60 


3.96 


13.75 


1 


.007 





4 


-1.21 


1.52 


4.57 


4.44 


14.49 


1 


.005 





3 


-0.75 


1.91 


5.32 


4.80 


15.08 


1 


.003 





2 


-0.37 


2.22 


5.92 


5.09 


15.60 


1 


.002 


0. 


1 


-0.11 


2.41 


6.33 


5.30 


15.85 


1 


.000 








0.00 


2.49 


6.47 


5.37 


15.99 


1 


.000 



TABLE IV: Values of reduced internal and excess free energies, chemical potential, pressure, and inverse compressibility /3 

as a function of the coverage x f° r a fixed state point p* — 0.68 and T* = 1.0. The last column gives the corresponding values 
of <To/a. 



X 


9HH (cr + ) 


5hh(Act ) 


3x(a+) 


,9x(Acr ) 


3ht(ct + ) 


#ht(Act ) 


z 


1.0 


2.476 


1.367 


2.476 


1.367 


2.476 


1.367 


10.64 


0.9 


2.574 


1.527 


2.499 


1.468 


2.449 


0.966 


9.05 


0.8 


2.745 


1.257 


2.631 


1.284 


2.781 


0.814 


7.53 


0.7 


3.287 


1.211 


3.092 


1.280 


2.367 


0.649 


6.04 


0.6 


4.101 


1.431 


3.906 


1.483 


2.052 


0.608 


4.75 


0.5 


4.643 


1.802 


4.907 


1.806 


1.941 


0.659 


3.53 


0.4 


4.402 


2.015 


2.094 


0.758 


1.986 


0.745 


2.42 


0.3 


4.025 


1.869 


2.219 


0.771 


2.093 


0.778 


1.51 


0.2 


4.429 


1.708 


2.341 


0.747 


2.298 


0.754 


0.74 


0.1 


6.199 


1.918 


2.678 


0.776 


2.668 


0.773 


0.22 


0.0 


3.069 


0.849 


3.069 


0.849 


3.069 


0.849 


0.0 



TABLE V: Behavior of g a {cr + ) and g a (\a~), where a = {HH,X,HT}. As in the text, HH, X, and HT refer to ni • n 2 = -1, 0, 1 
respectively. All values refer to the same state point p* = 0.68 and T* — 1.0. Here and below, a + and Act - refer to values at 
discontinuities inside the well, a < r < Act. The last column gives the average coordination number z within the well. 



X g(a + ,-l) ff(ACT~,-l) g(CT+,0) ff (Act~ , 0) g(CT+,+l) g(ACT-,+l) 



1.0 


2.476 


1.367 


2.476 


1.367 


2.476 


1.367 


0.9 


2.802 


1.343 


2.585 


1.264 


2.591 


1.263 


0.8 


3.073 


1.312 


2.650 


1.169 


2.653 


1.151 


0.7 


3.238 


1.289 


2.762 


1.101 


2.681 


1.032 


0.6 


3.373 


1.263 


2.869 


1.046 


2.558 


0.891 


0.5 


3.444 


1.220 


2.927 


0.989 


2.407 


0.775 


0.4 


3.513 


1.174 


2.969 


0.938 


2.513 


0.764 


0.3 


3.540 


1.114 


2.977 


0.891 


2.773 


0.810 


0.2 


3.480 


1.031 


2.955 


0.849 


2.926 


0.832 


0.1 


3.380 


0.952 


3.016 


0.842 


3.023 


0.844 


0.0 


3.069 


0.849 


3.069 


0.849 


3.069 


0.849 



TABLE VI: Behavior of g{cr + , cos 612) and <?(Act , cos #12) for cos #12 = —1,0, 1 corresponding to HH, X and HT configurations 
respectively. All values are as in the previous table. 
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Figure Captions 





HT X HH 



FIG. 1: The one-patch Kern-Frenkel potential, where fi2 is the direction joining the two centers (top panel). Directions of 
the patches are specified by unit vectors ni and n2 . Of particular relevance are the relative orientations of the two patches as 
defined by #12 = 82 ~ Oi' 812 = (HT), 812 — it/2 (X), and #12 = vr (HH) (bottom panel). Note that a larger, lighter spherical 
surface corresponds to the attractive part, meaning here a coverage \ > 0.5. 
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exp(Pn) 



FIG. 2: Plot of reduced density p* versus e"^ for \ — 80% coverage (top panel). Results are for both gas and liquid branches 
as computed from GCMC and RHNC. The bottom panel focuses only on the gas branch to enlarge the small-density region 
and display consistency with the correct ideal gas limit. Table HT1 provides the numerical values for both GEMC and RHNC 
calculations. 
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FIG. 3: Phase diagram for \ = 0.8 (80% coverage). Circled points indicate that an extrapolation procedure has been used in 
the RHNC data. Table Hill provides the RHNC numerical values. 
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T=1.00, p =0.10 



HT (8 I3 =0) 
i X (9 ]2 =ji/2) 
> HH (9 =71) 




1.5 



2 

r/a 



2.5 



FIG. 4: Plots of g(12) as a function of r/a for p* = 0.68 and temperatures T* = 0.65 (top panel) and T* = 1.0 (center panel). 
Different lines refer to three specific patch orientations having angles 812 = (HT), 612 = n/2 (X), and #12 = t (HH), all at 
X = 0.8 coverage. For comparison, in the bottom panel the corresponding low-density case p* = 0.1, T* = 1.0 is reported. SW 
and HS indicate the isotropic square-well and hard-sphere limits, respectively. Note that all panels have been drawn on the 
same scale. 
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FIG. 5: (color online) Snapshots of two 80% coverage configurations with p* — 0.68 and temperatures T* 
and T* = 1.0 (center panel). The last panel depicts a corresponding low-density case with p* = 0.1, T* 
selected state points are the same values considered in previous figures. 



= 0.65 (top panel) 
= 1.0. The three 



23 




FIG. 6: Plots of g(r, cos #12) as a function of r/a for p* = 0.68 and temperatures T* = 0.65 (top panel) and T* — 1.0 (center 
panel). All results pertain to \ = 0.8 coverage and both MC and RHNC results are reported. The bottom panel depicts the 
same quantities for the low-density case p* = 0.1, T* — 1.0. 
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FIG. 7: Plots of the coefficieunt Sooo(k) as a function of ka for p* = 0.68 and temperatures T* — 0.65 (top panel) and T* = 1.0 
(center panel). The bottom panel depicts the low-density case p* = 0.1, T* = 1.0. All results pertain to x — 0-8 coverage and 
both MC and RHNC results are reported. 
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FIG. 8: Plots of g(12) as a function oir/a for p* — 0.68 and temperature T* = 1.00 at different coverages ranging from \ — 1-0 
(fully symmetric square-well) to x = (hard spheres). Different panes refer to different configurations: HH (top panel), X 
(center panel), HT (bottom panel). Note that all curves have been drawn on the same scale. 
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FIG. 9: Plots of g(r,cos8) as a function of rja for p* = 0.68 and temperature T* = 1.00 at different coverages ranging from 
X = 1.0 (fully symmetric square-well) to x = (hard spheres). Different panels refer to different configurations as before: HH 
(top panel), X (center panel), HT (bottom panel). Again, all curves have been drawn on the same scale. 
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FIG. 10: Plots of Sboo(fc) as a function of ka for p* — 0.68 and temperature T* = 1.00 at different coverages ranging from 
X = 1-0 (fully symmetric square-well) to x — (hard spheres). 



